function [volume_calibration, use_average_image] = deform_images_psuedo_phase_toolbox(patient, phases, use_average_image)
    
    %| Deform the reference image to *phases* using the model parameters;
    %| The 5d motion model generates DVFs going from an arbitrary phase,
    %| to the reference phase, so these need to be inverted.
    
    %This code uses backwards2forwards.m (mex code compiled from C) from Mulitmodality non-rigid demon algorith imagef resgistration toolbox
    % written by Dirk-Jan Kroon (http://www.mathworks.com/matlabcentral/fileexchange/21451-multimodality-non-rigid-demon-algorithm-image-registration).
    
    %|  The DVFs are then written to file, and an Elastix
    %|  TransformParameters.0.txt file created (using the originial Elastix
    %|  registration as a template), so that the deformation can be
    %|  performed by Transformix.  There are many easy ways to do this
    %|  (3dinterp in Matlab for instance), but I found Transformix to be the
    %|  most accurate in deforming an image. - all this is done within the
    %|  'deform_using_transformix_toolbox.m' function.
    
    %|  Once 2 images at different static phases (default is 5th & 85th%ile), the
    %|  CT lung volume is calculated and a linear fit between volume and
    %|  bellows voltage calculated.  Lung volume from the psuedo CT images
    %|  are calculated from the lung mask (sum of all voxels), and also using
    %|  the partial are content of the image (more correct, but potentially more room for error).
    %|  The last line of this function chooses which to use.
    
    %|          *Please make sure Matlab can find your install of Elastix
    %|                  patient.elastix_root_folder needs to be set in
    %|                  choose_patient_toolbox.m
    %|                      default = "C:/Program Files/elastix/" -
    %|                                                          (windows).
    
    
    %------------------------------------------------------------------------
    %|      Dependancies;
    %|                  invert_motion_field
    %|                  deform_using_transformix_toolbox
    %|
    %------------------------------------------------------------------------
    %   This file is part of the
    %   5D-Novel4DCT Toolbox  ("Novel4DCT-Toolbox")
    %   DH Thomas, Ph.D
    %   University of California, Los Angeles
    %   Contact: mailto:dhthomas@mednet.ucla.edu
    %------------------------------------------------------------------------
    % $Author: DHThomas $	$Date: 2014/04/01 10:23:59 $	$Revision: 0.1 $
    
    
    if nargin < 3
        %         use_average_image = 0; %default - don't used the averaged image to deform (use the reference instead).
%         prompt = 'Would you like to deform the averaged image (1) or the reference image [0] (press return for 0)?';
        use_average_image = 0;%input(prompt);
    end
    
    if nargin < 2
        phases = [5,85];
    end
    
    
    
    fh=findall(0,'type','figure');
    fig_num = max(fh)+1;
    
    patient.folder_psuedo_image = [patient.model_folder '/Psuedo_images'];
    mkdir(patient.folder_psuedo_image)
    str = 'Y';
    if exist([patient.model_folder '/simple_coefs.mat'])>0
        
        prompt = 'Bellows Calibration exists. \n                Do you want to re-run Calibration code? Y [Hit Return for N]: ';
        str = input(prompt,'s');
        if isempty(str)
            str = 'N';
        end
    end
    if strcmpi('Y', str)>0;str = 'Y';end
    if str == 'N';
        load([patient.model_folder '/simple_coefs']);
        
    else
        display(sprintf('       Deforming reference image (%d) to phases; [%s]', patient.ref, num2str(phases)))
        
        if use_average_image>0;
            load([patient.folder '/img_average'])
            deform_image = img_average;
        else
            deform_image = patient.static;
        end
        
        load([patient.folder '/hu_fit'])
        hu_fit_p1=squeeze(hu_fit(1,:,:,:));
        hu_fit_p2=squeeze(hu_fit(2,:,:,:));
        clear hu_fit
        
        step=1;
        
        ii=1:patient.dim(1);
        jj=1:step:patient.dim(2);
        kk=1:step:patient.dim(3);
        
        %
        % bv_smooth=sort((patient.bellows_volt_drifted),'descend');
        % % inhale=bv_smooth(round(R85))
        % % exhale=bv_smooth(round(R15))
        % %
        
        volt_5pct=-prctile(-patient.bellows_volt_drifted(:),5);
        volt_85pct=-prctile(-patient.bellows_volt_drifted(:),85);
        
        [~,~,grid_z]=ndgrid(ii,jj,kk);
        volt_ref_img=reshape(patient.bellows_volt_drifted(grid_z,1),size(ii,2),size(jj,2),size(kk,2));
        
        display('               Loading model parameters...');
        load([patient.model_folder '/model_params/alpha']);
        load([patient.model_folder '/model_params/beta']);
        load([patient.model_folder '/model_params/constant']);
        model_u=single([ constant.u(:) alpha.u(:) beta.u(:)]);% gamma_vec_u];
        model_v=single([ constant.v(:) alpha.v(:) beta.v(:)]);% gamma_vec_u];
        model_w= single([ constant.w(:) alpha.w(:) beta.w(:)]);% gamma_vec_u];
        clear alpha beta constant
        
        for phase=1:size(phases,2)
            volt(phase)=-prctile(-patient.bellows_volt_drifted(:),phases(phase));
            flow(phase)=0;
            hu_fit_p2_finite=hu_fit_p2(isfinite(hu_fit_p2));
            hu_correction= nanmean(hu_fit_p2(patient.static_mask>0))*(volt(phase)-volt_ref_img);%ones(size(hu_fit_p2)).*
            %     hu_correction= (hu_fit_p2).*(volt(phase)-volt_ref_img);%ones(size(hu_fit_p2)).*
            
            figure(fig_num);hold on
            plot(patient.bellows_volt_drifted,'b--')
            plot([1,patient.dim(3)],[volt(phase),volt(phase)],'r--','LineWidth',2);
            hold off
            
            tic
            
            E=[1 volt(phase) flow(phase)];
            
            vX = double(reshape(model_u*E', size(ii,2), size(jj,2), size(kk,2)));
            vY = double(reshape(model_v*E', size(ii,2), size(jj,2), size(kk,2)));
            vZ = double(reshape(model_w*E', size(ii,2), size(jj,2), size(kk,2)));
            
            tic
            display(sprintf('               Deforming Image %d of %d', phase,size(phases,2)));
            [ivy,ivx,ivz]=backwards2forwards(vY,vX,vZ);
            ivy(ivy == 0) = NaN;
            ivx(ivx == 0) = NaN;
            ivz(ivz == 0) = NaN;
            image_ref = patient.static;%
            %         image_ref(patient.body_mask<1)= min(patient.static(:));
            image_scan_pseudo_tmp = deform_image_mirt3d(image_ref, ivy, ivx, ivz);
            image_scan_pseudo_tmp(isnan(image_scan_pseudo_tmp)) = min(image_scan_pseudo_tmp(:));
            mask_scan_pseudo = deform_image_mirt3d(patient.static_mask, ivy, ivx, ivz);
            mask_scan_pseudo(isnan(mask_scan_pseudo)) = 0;
            %         mask_scan_pseudo(isnan(mask_scan_pseudo)) = min(mask_scan_pseudo(:));
            deformed_mask = logical(mask_scan_pseudo);
            hu_correction(mask_scan_pseudo==0)=0;
            deformed_image=image_scan_pseudo_tmp+hu_correction;
            
            phase_folder = [patient.folder_psuedo_image sprintf('/phase%d',phases(phase))]
            
            clear deformed_mask_tmp;
            clear deformed_image_tmp;
            
            partial_air=(1-((deformed_image+943)/(44+943)));
            
            volume_simple(phase)=sum(deformed_mask(:))*patient.voxel_size*0.001/1000;
            volume_partial(phase)=sum(partial_air(deformed_mask>0))*patient.voxel_size*0.001/1000;
            
            
            
            save([patient.folder_psuedo_image sprintf('/Deformed_pseudo_Image_phase_%d',phases(phase))],'deformed_image')
            save([patient.folder_psuedo_image sprintf('/Deformed_pseudo_Mask_phase_%d',phases(phase))],'deformed_mask')
            
            if phases(phase)==5;
                
                image_X0=deformed_image;
                mask_X0=deformed_mask;
                save([patient.model_folder '/image_X0'], 'image_X0')
                save([patient.model_folder '/mask_X0'], 'mask_X0')
            end
            toc
            
        end
        
        
        
        simple_coefs= polyfit(volt,volume_simple,1)
        partial_coefs= polyfit(volt,volume_partial,1)
        
        figure;
        set(gcf,'units','normalized','position',patient.default_figure_position);
        
        hold on;
        p_s=plot(volt,volume_simple,'bo','LineWidth',2)
        f_s=plot(volt,polyval(simple_coefs,volt),'b--','LineWidth',2)
        p_p=plot(volt,volume_partial,'rd','LineWidth',2)
        f_p=plot(volt,polyval(partial_coefs,volt),'r-.','LineWidth',2)
        legend([f_s(1,1);f_p(1,1)],sprintf('Simple Volume \n Slope = %.2f L/V', simple_coefs(1)),sprintf('Partial Air Content \n Slope = %.2f L/V', partial_coefs(1)),'Location','E')
        xlabel('Bellows Voltage [V]');
        ylabel('Volume [L]')
        
        save([patient.model_folder '/volume_simple'], 'volume_simple');
        save([patient.model_folder '/volume_partial'], 'volume_partial');
        save([patient.model_folder '/simple_coefs'], 'simple_coefs');
        save([patient.model_folder '/partial_coefs'], 'partial_coefs');
        
        
    end
    
    volume_calibration = simple_coefs;
    %     volume_calibration = partial_coefs;
    
    
    
